
B.Est = function(data, cnt, K){
  
  N    = nrow(data)
  Lkp1 = as.vector(data[,(pL0 + K + 2) : (pL0 + 2*K + 2)])
  Akp1 = as.vector(data[,(pL0+1):(pL0+K+1)])
  # Ak   = cbind(rep(100,N),data[,(pL0+1):(pL0+K)])
  # Lk   = cbind(rep(100,N),data[,(pL0+K+2):(pL0+2*K+1)])
  # AkLk00 = as.numeric(Ak == 0 & Lk == 0)
  # AkLk01 = as.numeric(Ak == 0 & Lk == 1)
  # AkLk10 = as.numeric(Ak == 1 & Lk == 0)
  # AkLk11 = as.numeric(Ak == 1 & Lk == 1)
  # AkLk.  = as.numeric(Ak == 100 & Lk == 100)
  
  X.base = data[,1:pL0]
  temp   = list()
  for(i in 1:(K+1)) temp[[i]] = X.base
  X      = do.call(rbind, temp)
  
  X.model = model.matrix(~ 1 + Lkp1 + X )
  wt      = rep(cnt, K+1)
  
  b.model = glm.fit(X.model, Akp1, family=binomial(), weights = wt)
  
  return(b.model$coefficients)
  
}










